{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sisl\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analysing output from Siesta makes analysis much easier since many things are intrinsically enabled through `sisl` - i.e. orbital numbering etc. may easily be handled with `sisl`.\n",
    "\n",
    "In this example we will use Siesta to calculate the *heavy* part, i.e. the projected density-of-states for a fine $k$-point sampling. The formula used to calculate the projected density-of-states is:\n",
    "\\begin{align}\n",
    "  \\mathrm{PDOS}_\\nu(E) &= \\int \\mathrm d k \\sum_i \\psi^*_{i,k,\\nu} [\\mathbf S(k) | \\psi_{i,k}\\rangle]_\\nu D(E-\\epsilon_{i}(k)),\n",
    "  \\\\\n",
    "  \\mathrm{DOS}(E) &= \\int \\mathrm d k \\sum_i D(E-\\epsilon_{i}(k)) = \\sum_\\nu \\mathrm{PDOS}_\\nu(E),\n",
    "\\end{align}\n",
    "Where *E* is the energy, $D(\\Delta E)$ is the distribution function (in Siesta, this is the Gaussian distribution function) and $\\epsilon_{i}(k)$ and $\\mathbf S(k)$ are the ith eigenvalue and overlap matrix, respectively, at a given **k**-point. $\\psi_{i,k}$ is the ith **k**-resolved eigenstate, and $\\psi^*_{i,k,\\nu}$ is the corresponding complex-conjugate resolved for a concrete atomic orbital $\\nu$. \n",
    "Finally, we will do the same calculation using `sisl` to ensure the correctness of both methods.\n",
    "Note that `sisl` is (for now) not MPI parallelized and thus for large systems it may be *way* more efficient to use Siesta to calculate PDOS, and then use `sisl` to post-process the PDOS data.\n",
    "\n",
    "The system will again be graphene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene = sisl.geom.graphene(1.44)\n",
    "graphene.write(\"STRUCT.fdf\")\n",
    "graphene.write(\"STRUCT.xyz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To tell Siesta to calculate the projected DOS, one has to add these flags to the input (if you want you can play with the numbers):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "open(\"PDOS.fdf\", \"w\").write(\n",
    "    \"\"\"\n",
    "# K-point sampling\n",
    "%block PDOS.kgrid.MonkhorstPack\n",
    "  63  0 0\n",
    "   0 63 0\n",
    "   0  0 1\n",
    "%endblock\n",
    "# Energy grid to calculate the PDOS on\n",
    "# E_min = -20 eV\n",
    "# E_max = 15 eV\n",
    "# Gaussian width: \\sigma 2 ^{1/2} = 0.2 eV\n",
    "# Number of points: 3500 (spacing of 0.01 eV)\n",
    "%block Projected.DensityOfStates\n",
    "-20.00 15.00 0.200 3500 eV\n",
    "%endblock\n",
    "\"\"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please open `PDOS.fdf` and `RUN.fdf`.\n",
    "- Why do you see two different $k$-point samplings?\n",
    "\n",
    "Now run Siesta to calculate the electronic structure and the PDOS:\n",
    "\n",
    "     siesta RUN.fdf | tee RUN.out\n",
    "\n",
    "Subsequently, we will read in the PDOS information and post-process it.  \n",
    "First we will extract data found in the `siesta.PDOS.xml` file. Three quantities will be returned:\n",
    "1. The geometry (as found in the XML file) - this also contains all atomic species etc.\n",
    "2. The energy grid used\n",
    "3. Orbital and energy resolved DOS\n",
    "\n",
    "Before analysing, please examine the data returned by viewing the shapes and printing the geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geometry, E, PDOS = sisl.get_sile(\"siesta.PDOS.xml\").read_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Examine output by printing shapes and geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "1. Plot the total DOS (invoke Eq. 2 in the top paragraph)\n",
    "2. Plot the total PDOS on the first atom  \n",
    "   *HINT*: Sum all orbitals on the first atom\n",
    "3. Figure out which orbital contributes to the Dirac cone by plotting the PDOS for all orbitals\n",
    "4. Read in the Hamiltonian (as done in [S 1](../S_01/run.ipynb)) and create a Monkhorst-Pack grid as the one in Siesta input. Use this snippet to calculate the PDOS using `sisl`\n",
    "\n",
    "    ```python\n",
    "    mp = sisl.MonkhorstPack(H, [nx, ny, 1])\n",
    "    def wrap_multiple(eigenstate):\n",
    "        DOS = eigenstate.DOS(E)\n",
    "        PDOS = eigenstate.PDOS(E)\n",
    "        return sisl.oplist([DOS, PDOS])\n",
    "    DOS, PDOS = mp.apply.average.eigenstate(wrap=wrap_multiple, eta=True)\n",
    "    \n",
    "    ```\n",
    "\n",
    "    Search the API documentation for the `MonkhorstPack.apply.average` method and figure out what it does. Note, there are other `MonkhorstPack.apply.*` methods - these are all extremely useful when calculating many quantities of data in a Brillouin-zone object.\n",
    "5. Ensure that the `DOS` and the summed `PDOS` (use Eq. 2) are equivalent.\n",
    "6. Compare the Siesta PDOS with sisl PDOS, why are they different?  \n",
    "   *HINT*: $\\sigma$. Check the API for the `PDOS` method in the Siesta manual, provided via this tutorial (`../siesta.pdf`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(E, PDOS.sum(0), label=\"DOS\")\n",
    "plt.xlabel(r\"$E - E_F$ [eV]\")\n",
    "plt.ylabel(r\"DOS [1/eV]\")\n",
    "\n",
    "# Add more plots for individual PDOS\n",
    "\n",
    "\n",
    "# Add legend(s)\n",
    "plt.legend();"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
